library(adehabitat)
library(fields)

#create a pyramid
tnum = 12 #the width of the base of the square pyramid
pdata = matrix(NA,nrow=tnum,ncol=tnum)
for (i in 1:tnum) {
  for (j in 1:tnum){
    if (i <= tnum/2 & j <= tnum/2) {
      if (j < i) {
        pdata[i,j] = j
        pdata[tnum-i+1,j] =  j
        pdata[tnum-i+1,tnum-j+1] =  j
        pdata[i,tnum-j+1] = j
      } else {
        pdata[i,j] = i
        pdata[tnum-i+1,j] =  i
        pdata[tnum-i+1,tnum-j+1] =  i
        pdata[i,tnum-j+1] =  i
      }
    }
  }
}
pdata
pdata = rbind(pdata,pdata)
################################################################################################
png("c:/R/course/data/mde.png",width = 2000,height=1200)
par(mfrow=c(3,3),mar=c(0,1,2,1),cex=1,cex.lab=1.5,cex.main=1.75)
Latitude = c(1:(tnum*2))
Longitude = c(1:tnum)
Elevation = pdata
################################################################################################
#create the 1D longitude
#z2 given x between 3 and 8
z2 = pdata
for (i in 1:tnum){
  if (i >=3 & i<=8) {
    z2[,i] = 1
  } else {
    z2[,i] = 0
  }
}
#create the plot and bounding box
persp(x = Latitude,y = Longitude, z=Elevation, scale=F, box=T, theta=45, phi = 20, col=drape.color(z2,col=c("grey70","grey30")), border="white")->pmat
title(main = "Longitudinal Range")
p1 = trans3d(24,3,1,pmat)
p2 = trans3d(24,8,1,pmat)
p3 = trans3d(24,8,6,pmat)
p4 = trans3d(24,3,6,pmat)
p5 = trans3d(1,3,1,pmat)
p6 = trans3d(1,8,1,pmat)
p7 = trans3d(1,8,6,pmat)
p8 = trans3d(1,3,6,pmat)
lines(rbind(p1,p2,p3,p4,p1,p5,p6,p7,p8,p5),lty=2)
lines(rbind(p4,p8),lty=2)
lines(rbind(p3,p7),lty=2)
lines(rbind(p2,p6),lty=2)
################################################################################################
#create the 1D latitude
#z2 given x between 4 and 9
z2 = pdata
for (i in 1:(tnum*2)){
  if (i >=5 & i<=21) {
    z2[i,] = 1
  } else {
    z2[i,] = 0
  }
}
#create the plot and bounding box
persp(x = Latitude,y = Longitude, z=Elevation, scale=F, box=T, theta=45, phi = 20, col=drape.color(z2,col=c("grey70","grey30")), border="white")->pmat
title(main = "Latitudinal Range")
p1 = trans3d(21,1,1,pmat)
p2 = trans3d(21,12,1,pmat)
p3 = trans3d(21,12,6,pmat)
p4 = trans3d(21,1,6,pmat)
p5 = trans3d(5,1,1,pmat)
p6 = trans3d(5,12,1,pmat)
p7 = trans3d(5,12,6,pmat)
p8 = trans3d(5,1,6,pmat)
lines(rbind(p1,p2,p3,p4,p1,p5,p6,p7,p8,p5),lty=2)
lines(rbind(p4,p8),lty=2)
lines(rbind(p3,p7),lty=2)
lines(rbind(p2,p6),lty=2)
################################################################################################
#create the 1D elevation
#z2 given x between 2 and 4
z2 = pdata
for (i in 1:tnum){
  for (j in 1:(tnum*2)) {
    if (pdata[j,i] >=2 & pdata[j,i]<=5) {
      z2[j,i] = 1
    } else {
      z2[j,i] = 0
    }
  }
}
z2[8,5] = 0;z2[8,8] = 0;z2[5,5] = 0;z2[5,8] = 0;z2[17,5] = 0;z2[17,8] = 0;z2[20,5] = 0;z2[20,8] = 0
#create the plot and bounding box
persp(x = Latitude,y = Longitude, z=Elevation, scale=F, box=T, theta=45, phi = 20, col=drape.color(z2,col=c("grey70","grey30")), border="white")->pmat
title(main = "Elevational Range")
p1 = trans3d(24,1,2,pmat)
p2 = trans3d(24,12,2,pmat)
p3 = trans3d(24,12,5,pmat)
p4 = trans3d(24,1,5,pmat)
p5 = trans3d(1,1,2,pmat)
p6 = trans3d(1,12,2,pmat)
p7 = trans3d(1,12,5,pmat)
p8 = trans3d(1,1,5,pmat)
lines(rbind(p1,p2,p3,p4,p1,p5,p6,p7,p8,p5),lty=5)
lines(rbind(p4,p8),lty=5)
lines(rbind(p3,p7),lty=5)
lines(rbind(p2,p6),lty=5)
################################################################################################
#create the 2D longitude x latitude
#z2 given x between 2 and 4
z2 = pdata
for (i in 1:tnum){
  for (j in 1:(tnum*2)) {
    if (i >=3 & i<=8 & j >=5 & j<=21) {
      z2[j,i] = 1
    } else {
      z2[j,i] = 0
    }
  }
}

#create the plot and bounding box
persp(x = Latitude,y = Longitude, z=Elevation, scale=F, box=T, theta=45, phi = 20, col=drape.color(z2,col=c("grey70","grey30")), border="white")->pmat
title(main = "Long x Lat Range")
p1 = trans3d(21,3,1,pmat)
p2 = trans3d(21,8,1,pmat)
p3 = trans3d(21,8,6,pmat)
p4 = trans3d(21,3,6,pmat)
p5 = trans3d(5,3,1,pmat)
p6 = trans3d(5,8,1,pmat)
p7 = trans3d(5,8,6,pmat)
p8 = trans3d(5,3,6,pmat)
lines(rbind(p1,p2,p3,p4,p1,p5,p6,p7,p8,p5),lty=5)
lines(rbind(p4,p8),lty=5)
lines(rbind(p3,p7),lty=5)
lines(rbind(p2,p6),lty=5)
################################################################################################
#create the 2D longitude x elevation
#z2 given x between 2 and 4
z2 = pdata
for (i in 1:tnum){
  for (j in 1:(tnum*2)) {
    if (i >=3 & i<=8 & pdata[j,i] >=2 & pdata[j,i]<=5) {
      z2[j,i] = 1
    } else {
      z2[j,i] = 0
    }
  }
}
z2[8,5] = 0;z2[8,8] = 0;z2[5,5] = 0;z2[5,8] = 0;z2[17,5] = 0;z2[17,8] = 0;z2[20,5] = 0;z2[20,8] = 0
#create the plot and bounding box
persp(x = Latitude,y = Longitude, z=Elevation, scale=F, box=T, theta=45, phi = 20, col=drape.color(z2,col=c("grey70","grey30")), border="white")->pmat
title(main = "Long x Elev Range")
p1 = trans3d(24,3,2,pmat)
p2 = trans3d(24,8,2,pmat)
p3 = trans3d(24,8,5,pmat)
p4 = trans3d(24,3,5,pmat)
p5 = trans3d(1,3,2,pmat)
p6 = trans3d(1,8,2,pmat)
p7 = trans3d(1,8,5,pmat)
p8 = trans3d(1,3,5,pmat)
lines(rbind(p1,p2,p3,p4,p1,p5,p6,p7,p8,p5),lty=5)
lines(rbind(p4,p8),lty=5)
lines(rbind(p3,p7),lty=5)
lines(rbind(p2,p6),lty=5)
################################################################################################
#create the 2D latitude x elevation
#z2 given x between 2 and 4
z2 = pdata
for (i in 1:tnum){
  for (j in 1:(tnum*2)) {
    if (j >=5 & j<=21 & pdata[j,i] >=2 & pdata[j,i]<=5) {
      z2[j,i] = 1
    } else {
      z2[j,i] = 0
    }
  }
}
z2[8,5] = 0;z2[8,8] = 0;z2[5,5] = 0;z2[5,8] = 0;z2[17,5] = 0;z2[17,8] = 0;z2[20,5] = 0;z2[20,8] = 0
#create the plot and bounding box
persp(x = Latitude,y = Longitude, z=Elevation, scale=F, box=T, theta=45, phi = 20, col=drape.color(z2,col=c("grey70","grey30")), border="white")->pmat
title(main = "Lat x Elev Range")
p1 = trans3d(21,1,2,pmat)
p2 = trans3d(21,12,2,pmat)
p3 = trans3d(21,12,5,pmat)
p4 = trans3d(21,1,5,pmat)
p5 = trans3d(5,1,2,pmat)
p6 = trans3d(5,12,2,pmat)
p7 = trans3d(5,12,5,pmat)
p8 = trans3d(5,1,5,pmat)
lines(rbind(p1,p2,p3,p4,p1,p5,p6,p7,p8,p5),lty=5)
lines(rbind(p4,p8),lty=5)
lines(rbind(p3,p7),lty=5)
lines(rbind(p2,p6),lty=5)
################################################################################################
#create the 3D longitudexlatitude x elevation
#z2 given x between 2 and 4
z2 = pdata
for (i in 1:tnum){
  for (j in 1:(tnum*2)) {
    if (i >=3 & i<=8 & j >=5 & j<=21 & pdata[j,i] >=2 & pdata[j,i]<=5) {
      z2[j,i] = 1
    } else {
      z2[j,i] = 0
    }
  }
}
z2[8,5] = 0;z2[8,8] = 0;z2[5,5] = 0;z2[5,8] = 0;z2[17,5] = 0;z2[17,8] = 0;z2[20,5] = 0;z2[20,8] = 0
#create the plot and bounding box
persp(x = Latitude,y = Longitude, z=Elevation, scale=F, box=T, theta=45, phi = 20, col=drape.color(z2,col=c("grey70","grey30")), border="white")->pmat
title(main = "Long x Lat x Elev Range")
p1 = trans3d(21,3,2,pmat)
p2 = trans3d(21,8,2,pmat)
p3 = trans3d(21,8,5,pmat)
p4 = trans3d(21,3,5,pmat)
p5 = trans3d(5,3,2,pmat)
p6 = trans3d(5,8,2,pmat)
p7 = trans3d(5,8,5,pmat)
p8 = trans3d(5,3,5,pmat)
lines(rbind(p1,p2,p3,p4,p1,p5,p6,p7,p8,p5),lty=5)
lines(rbind(p4,p8),lty=5)
lines(rbind(p3,p7),lty=5)
lines(rbind(p2,p6),lty=5)
################################################################################################
dev.off()
################################################################################################
